Steven Dutch, Professor Emeritus, Natural and Applied Sciences, Universityof Wisconsin - Green Bay
7-1 DIP AND STRIKE
7-1a ATTITUDE COSINES OF A DIPPING PLANE
Imagine a line of unit length perpendicular to a horizontalplane. The plane lies in the Z = 0 plane and the line (the poleto the plane) points along the positive Z-axis. Now rotate theplane around a north-south axis through a dip angle D. Obviouslythe pole to the plane rotates also. Its projection along theZ-axis is COS(D) and its projection in the Z = 0 plane is SIN(D).
Right away we encounter an ambiguity. The plane can have oneof two dip directions. We can remove this ambiguity by definingdip according to a convention. If we perform our dip rotationtoward the east, the projection of the pole lies along the+X-axis, which is convenient for computation. The strike can beeither 0 or 180. Let us define the strike in such a way that whenwe look along in the strike direction, the plane dips down and toour right. Thus a plane of strike 0 dips east and one of strike180 dips west. It is apparent that if the strike is S, thedip direction as we have defined it is S + 90.
Now rotate the plane around the Z-axis by the strike angleS. The projection of the pole along the Z-axis is still COS(D).Its projections on the other axes require a bit of thought,however, because we conventionally place north at the top of mostmaps and measure direction clockwise. On the other hand, we putthe X-axis of Cartesian coordinates horizontally and measuredirections in a counter-clockwise sense. Thus we can relate thetwo coordinate systems as follows:
Compass
Cartesian Polar Coordinates
North 0/360
+Y axis
90
East 90
+X axis
0/360
South 180
-Y axis
270
West 270
-X axis
180Recall also that the projection of the pole in the Z = 0 plane,according to the convention we defined above, has the compassdirection S + 90. The projection of the pole on the X-axis willtherefore be SIN(D)*COS(S), and the projection along the Y-axiswill be -SIN(D)*COS(S). Therefore, the direction cosines of aplane with dip D and strike S are:
SIN(D)*COS(S), -SIN(D)*SIN(S), COS(D)It is easy to verify that the squares of the direction cosinessum to 1.
In order to help visualize these results, it is useful toconsider a few examples. For planes which strike along thecardinal directions we have:
Strike
Dip
Direction Cosines
NORTH
EAST
SIN(D)
0
COS(D)
EAST
SOUTH
0
-SIN(D) COS(D)
SOUTH
WEST
-SIN(D)
0
COS(D)
WEST
NORTH
0
SIN(D) COS(D)For oblique strikes, we see from the above that
STRIKE QUADRANT DIP QUADRANT SIGNS OF DIRECTION COSINES
NE
SE
+ - +
SE
SW
- - +
SW
NW
- + +
NW
NE
+ + +
If this sounds complicated, consider the computing problemsinvolved in a more traditional approach to strike and dip. Let'ssay, for example, you have written a program to accept thefollowing input:
INPUT STRIKE ? N 35 W
INPUT DIP ? 45 NEFirst of all, you have to write the program so that it caninterpret the letters and numbers. Then, and this is the complexpart, you have to provide for all the possible dip cases. The dipcould just as plausibly have been written 45 N or 45 E; theopposite dip direction could be written 45 SW, 45 S or 45 W. Andwhat if, as sometimes happens, you get your field notes garbledand input the dip as 45 NW or 45 SE ?7-1b A PROGRAM FOR READING TRADITIONAL ATTITUDE FORMATS
For those who prefer to be able to use the traditionalformat, here is a program to turn strike and dip into directioncosines.
10 C1=180/3.14159265359
100 INPUT"INPUT STRIKE ";A$,S,B$
101 IF A$="N" THEN 106
102 IF A$="S" THEN 106
103 GOTO 108
106 IF B$="E" THEN 110
107 IF A$="W" THEN 110
108 PRINT "THAT FORMAT NOT ALLOWED"
109 GOTO 100
110 INPUT"INPUT DIP ";D,C$
120 REM CASE 1--A$ = SOUTH
130 IF A$="N" THEN 200
140 A$="N":S=S+180
150 IF B$="E" THEN 180
160 B$="E"
170 GOTO 200
180 B$="W"
190 REM CHANGE TO CONVENTIONAL N-STRIKE-E/W FORMAT
200 REM A$ = NORTH -- CONVENTIONAL CASE
205 IF B$="E" THEN 220
210 S=360-S
220 REM SPECIAL CASES-CARDINAL DIRECTIONS
225 IF ABS(SIN(2*S/c1))> 0.000001 THEN 350
230 REM NORTH-SOUTH STRIKE
235 IF ABS(SIN(S/c1))>.001 THEN 300
245 IF C$="E" THEN 260
250 IF C$="W" THEN 260
255 GOTO 1030
260 IF C$="W" THEN 280
265 REM EAST DIP
270 S = 0
275 GOTO 1100
280 REM WEST DIP
285 S=180
290 GOTO 1100
300 REM EAST-WEST STRIKE
305 IF ABS(COS(S/c1))>.001 THEN 350
310 IF C$="N" THEN 325
315 IF C$="S" THEN 325
320 GOTO 1030
325 IF C$="S" THEN 335
328 REM NORTH DIP
330 S = 270
332 GOTO 1100
335 REM SOUTH DIP
340 S= 90
345 GOTO 1100
350 REM INTERPRET DIP DIRECTIONS
355 REM CASE 1-STRIKE INTO FIRST QUADRANT
360 IF COS(S)<0 THEN 500
365 IF SIN(S)<0 THEN 500
370 IF C$="NE" THEN 1030: REM ERROR
375 IF C$="SW" THEN 1030: REM ERROR
380 IF C$="S" THEN 410
385 IF C$="E" THEN 410
390 IF C$="SE" THEN 410
395 REM DIP IS NORTHWEST-STRIKE CONVENTION NOT FOLLOWED
400 S=S+180
410 REM STRIKE CONVENTION FOLLOWED-NO CHANGE
420 GOTO 1080
500 REM CASE 2-STRIKE INTO SECOND QUADRANT
520 IF COS(S)>0 THEN 700
530 IF SIN(S)<0 THEN 700
540 IF C$="SE" THEN 1030: REM ERROR
550 IF C$="NW" THEN 1030: REM ERROR
560 IF C$="S" THEN 610
570 IF C$="W" THEN 610
580 IF C$="SW" THEN 610
590 REM DIP IS NORTHEAST-STRIKE CONVENTION NOT FOLLOWED
600 S=S+180
610 REM STRIKE CONVENTION FOLLOWED-NO CHANGE
620 GOTO 1080
700 REM CASE 3-STRIKE INTO THIRD QUADRANT
720 IF COS(S)>1 THEN 900
730 IF SIN(S)>1 THEN 900
740 IF C$="NE" THEN 1030: REM ERROR
750 IF C$="SW" THEN 1030: REM ERROR
760 IF C$="N" THEN 810
770 IF C$="W" THEN 810
780 IF C$="NW" THEN 810
790 REM DIP IS SOUTHEAST-STRIKE CONVENTION NOT FOLLOWED
800 S=S+180
810 REM STRIKE CONVENTION FOLLOWED-NO CHANGE
820 GOTO 1080
900 REM CASE 4-STRIKE INTO FOURTH QUADRANT
920 IF COS(S)<1 THEN 1050
930 IF SIN(S)>1 THEN 1050
940 IF C$="NW" THEN 1030: REM ERROR
950 IF C$="SE" THEN 1030: REM ERROR
960 IF C$="N" THEN 1010
970 IF C$="E" THEN 1010
980 IF C$="NE" THEN 1010
990 REM DIP IS SOUTHWEST-STRIKE CONVENTION NOT FOLLOWED
1000 S=S+180
1010 REM STRIKE CONVENTION FOLLOWED-NO CHANGE
1020 GOTO 1080
1030 REM ERROR MESSAGE
1040 PRINT
1050 PRINT "DIP DIRECTION INCOMPATIBLE WITH STRIKE"
1060 PRINT
1070 GOTO 100
1080 REM COMPUTE DIRECTION COSINES
1100 C1 = SIN(D)*COS(S)
1110 C2 = -SIN(D)*SIN(S)
1120 C3 = COS(D)
1130 REM CONTINUE WITH MAIN PROGRAMOnce you have this subroutine working, it is a relatively simplematter to use it in any programs that involve inputtingstructural data. Just make a copy of an existing program thatemploys it; delete whatever you don't need, and begin building onthe existing foundation.
This program is a good example of how simple mathematicsgets built up into very long programs by IF statements. Theprogram is basically nothing but string comparison and somesimple spatial reasoning. A few comments are in order. Theconstant C1 defined in line 10 is necessary whenever you havedegrees and radians in the same program. The test in line 225 isan example of testing several things at once. Rather than test tosee if SIN(X) =0 for north-south strike and then test to see ifCOS(X)=0 for east-west strike, we can find one function thattests both conditions at the same time. The ABS function thatappears in several places is necessary because we cannot beentirely sure that, say, SIN(180/c1) will come out to be exactly equal to zero. We did, after all, truncate pi in theexpression for C1 and the computer may or may not yield a verysmall non-zero value.
The repetitious structure of the program lends itself nicelyto use of a good editor, so that blocks of program need bewritten only once, then copied and modified as desired. It iseasy to see, for example, that lines 350-420, 500-620, 700-820and 900-1020 are clones of one another apart from smallvariations. Be wary, though, of letting this sort of programmingget so mechanical that errors slip by.7-1C VARIATIONS IN USAGE
There are a few variations in terminology that need to benoted. First, we sometimes encounter the old term hade , whichis dip measured from the vertical rather than the horizontal, andis obviously (90 - dip). British maps sometimes describe planarstructures in terms of the trend and plunge of the down-dipdirection; that is a bed "dips 34 degrees toward 234". Theconversion to the convention we defined earlier isstraightforward. Simply subtract 90 degrees from the dipdirection to get the strike.7-2 TREND AND PLUNGE7-2a THE ORDINARY CASE
Trend and plunge are a lot easier to deal with because thereis no ambiguity in most cases. The plunge is the angle a linemakes with the horizontal. If the plunge is P and we consider theline to be plunging downward, the projection along the negativeZ-axis is obviously -SIN(P). The trend (call it T) is the strikeof a vertical line that contains the plunging line, and isspecified to be pointing in the direction of the plunge. In thiscase N20E is not the same as S20W, and 020 not the same as 200.
Suppose for the moment we specify trend in terms of azimuthfrom 0 to 360 degrees. The projection of the plunging line on ahorizontal plane is proportional to COS(P). The projection of theplunging line on the X-axis, then, is proportional toCOS(P)*SIN(T), and on the Y-axis to COS(P)*COS(T). Therefore, thedirection cosines of a line with trend T and plunge P are:
COS(P)*SIN(T), COS(P)*COS(T), -SIN(P)7-2b DEALING WITH UPWARD-DIRECTED LINES
In some cases, we must deal with directed lines in which itmakes a difference whether or not we consider the line to beplunging down or up. A measurement of sole-markings on a tiltedbed for paleo-current studies is a good example. Since the thirddirection cosine is -SIN(P), the simplest approach by far is torecord trend as usual, and define downward plunge (the usualcase) as positive, with upward plunge (only rarely needed) asnegative.7-2C A PROGRAM FOR READING TRADITIONAL TREND AND PLUNGE FORMATS
Just as in the case of strike and dip, it is useful to havea routine for reading quadrant-format field data.
10 C1=180/3.14159265359
100 INPUT"INPUT TREND ";A$,T,B$
101 IF A$="N" THEN 106
102 IF A$="S" THEN 106
103 GOTO 108
106 IF B$="E" THEN 110
107 IF A$="W" THEN 110
108 PRINT "THAT FORMAT NOT ALLOWED"
109 GOTO 100
110 INPUT"INPUT PLUNGE ";P
120 REM DEAL WITH EAST-WEST SENSE
130 IF B$="E" THEN 150
140 T=-T
150 DEAL WITH NORTH-SOUTH CONVENTION
160 IF A$="N" THEN 200
170 T=180-T
200 REM CALCULATE DIRECTION COSINES
210 K1= COS(P/c1)*SIN(T/c1)
220 K2= COS(P/c1)*COS(T/c1)
230 K3= -SIN(P/c1)Note how much simpler this program is than the correspondingprogram for strike and dip. As written, there is no need for aspecial provision for upward-directed lines; they can be treatedas having negative plunge.7-3 PITCH
It often happens that a linear structure will occur within aplane, such as a ripple mark on a bed or a mineral lineation on afoliation. It may be more convenient to describe the orientationof the line in terms of its orientation on the dipping planerather than in terms of trend and plunge, or we may want to refera line of known trend and plunge to a coordinate system based onthe dipping plane. In cases of this sort, we refer to the pitch of the line, which is the angle the line makes with the strikeline of the dipping plane as measured in that plane.
There is an immediate source of ambiguity involved. Say abed strikes north-south and dips east. A line with a pitch of 30degrees could plunge northeast if pitch is measured from thenorth end of the strike line, or southeast if pitch is measuredfrom the south end of the strike line.
We need some sort of convention. Since we're dealing withcompass directions in the field, define pitch angle in theclockwise (compass) sense. The pitch of any line that plungesstraight down-dip is always 90 degrees. In our example above,then, the northeast-plunging line has a pitch of 30 degrees andthe southeast-plunging line a pitch of 150 degrees. It is alsoapparent that a line of pitch zero will be pointing in the strikedirection as we defined it in 7-1a, so that our conventions aremutually consistent.
Developing a formula for pitch is a problem that is hard ifwe attempt a brute-force solution but very simple if we adopt abit of subtlety. In this case, project the dipping plane and theline onto a sphere. The great circle formed by the dipping planemakes an angle D with the equator of the sphere, where D, ofcourse, is the dip. Angle P1, the pitch,is the distance measured along thedipping great circle to the line. If wedraw a vertical great circle through theline, the distance P from the equatorto the dipping great circle is simply the plunge of the line, andthe distance T1 along the equator to the meridian is useful forfinding the trend of the line: T = S + T1. We have a rightspherical triangle, whose solution is elementary.
SIN(P) = SIN(P1)*SIN(D)
SIN(T1) = TAN(P)/TAN(D)
T = S + T17-4 ATTITUDES FROM DIRECTION COSINES
We frequently will need to recover the attitudes of a lineor plane from its direction or attitude cosines. Direction orattitude numbers can always be converted into the cosine form bynormalizing the sum of their squares to 1. Thus if we have aplane whose direction numbers are A, B, and C, the correspondingdirection cosines are A/SQR(A*A + B*B + C*C), and so on. Withthis point established, we can assume from this point on that weare dealing with direction cosines.7-4a STRIKE AND DIP FROM ATTITUDE COSINES
The attitude cosines of a plane are related to itsstrike S and dip D by the formula
A = SIN(D)*SIN(S), B = SIN(D)*COS(S), C = COS(D)If we have the attitude cosines A, B, and C, we have to workbackward. Obviously, D = K * ARCCOS(C), where K is ourradian-degree conversion factor. Also TAN(S) = A/b, but itdoes not follow that S = ATN(A/b), because S could also bedifferent by 180 degrees. The ideal function for convertingCartesian to Polar coordinates is ARCCOS:
S = K * (ARCCOS(B/SIN(D))*SGN(A)The principal value of ARCCOS is defined between 0 and 180. TheSGN function allows us to cover the range -180 to 0. If you wishto express angles in the range 0-360, the extra statements belowwill do the job.
100 IF S >= 0 THEN 120
110 S = S + 180
120 REM END EXAMPLE
If you are not fortunate enough to have a the ARCCOSfunction available, you must make use of ATN instead. The dip canbe found by D = K * ATN( SQR(1-C*C) / C). We can make use of thefact that TAN(S) = A/b if we can find some way to eliminate theambiguity between first and third quadrant solutions. Thefunction ATN is defined on the interval -90 to +90 degrees. Inthis interval COS(S) is always positive. If COS(S) is negative,then S differs from ATN(A/b) by 180 degrees. The followingformula will give a unique value of S
S = K * ATN(A/b) + 90 * (1-SGN(B))*SGN(A)The resulting value of S will lie in the range -180 to 180.If A or B are exactly zero the formula yields an erroneousresult, so any program that utilizes this formula must safeguardagainst such an error.7-4b TREND AND PLUNGE FROM DIRECTION COSINES
The formulas that relate trend and plunge to directioncosines are very similar in form to those for strike and dip, sothe solution is little more than a modification of the precedingcase. A line of trend T and plunge P has the direction cosines
A = SIN(T)*COS(P), B = COS(T)*COS(P), C = - SIN(P)If you have the full range of inverse trigonometric functions atyour disposal, the formulas are
P = - K * ARCSIN(C), T = K * (ARCCOS(B/cOS(P))*SGN(A)If you are restricted to the ATN function, the formulas become
P = K * ATN (C / SQR(1-C*C))
T = K * ATN(A/b) + 90 * (1-SGN(B))*SGN(A)Again K is our radian-degree conversion factor. Angles arecomputed in the range -180 to +180 and can be converted to therange 0-360 by the method shown in 7-47-5 CONVERTING NUMERICAL DATA INTO CONVENTIONAL FORMAT
Dip and plunge figures can be used just the way the programdelivers them. Strike and trend figures can as well if it ispermissible to give the figures in terms of azimuth from zero to360 degrees, and if the dip-direction convention is followed. Ifthe results are desired in the traditional quadrant format, oryou want to specify dip direction, a bit more work is needed.7-5a. TREND OR STRIKE
If you have access to the ARCCOS function, the followingsteps will work
10 K = 180/3.14159265359
......
1000 T = K * ARCCOS(B)
1100 T$ = "N " + STR$(ABS(T))
1200 IF A<0 THEN 1500
1300 T$ = T$ + " E"
1400 GOTO 1600
1500 T$ = T$ + " W"
1600 REM END EXAMPLEIn this case, we used T for trend but exactly the same procedurewill work for strike. If you are restricted to ATN, this set ofsteps will work for either strike or trend.
10 K = 180/3.14159265359
......
1000 T = K * ATN(A/b) - 90 * (1-SGN(A))
1100 T$ = "N " + STR$(ABS(T))
1200 IF A<0 THEN 1500
1300 T$ = T$ + " E"
1400 GOTO 1600
1500 T$ = T$ + " W"
1600 REM END EXAMPLE7-5b DIP DIRECTION
Since our direction cosines are defined in terms of ourstrike convention (7-1a), the dip direction is 90+S. Assuming wehave gotten numerical values for the strike and dip, we canexpress the dip direction in the traditional format as follows.In the example, S is the strike, D the dip ,D0 the dip direction(all in radians) and D$ the output string for the dip result.
10 K=180/3.1415926535
.......
1000 D0= S+90/c1
1010 D$="DIP = " + STR$(S*C1) + " "
1020 REM TEST DIP DIRECTIONS
1030 IF COS(S)<0 THEN 1070
1040 REM NORTH DIP COMPONENT
1050 D$ = D$ + "N"
1060 GOTO 1100
1070 REM SOUTH DIP COMPONENT
1080 D$ = D$ + "S"
1100 IF SIN(S)<0 THEN 1140
1110 REM EAST DIP COMPONENT
1120 D$ = D$ + "E"
1130 GOTO 1200
1140 REM WEST DIP COMPONENT
1150 D$ = D$ + "W"
1200 REM END EXAMPLE7-5 BEST FIT OF A PLANE (GREAT CIRCLE) TO A SET OF DIRECTIONS
This problem has many applications in geology. Thebest-known example is the fit of a fold axis to a collection ofbedding-plane attitudes. Another application is finding the poleof rotation for the opening of an ocean basin, given matchingpoints on opposite sides. We have a set of direction cosines,C1, C2, and C3 for each line. We'll assume the data is beinginput by the user, but READing the data from a set of DATAstatements is an easy variation. We want to find either a planethat is as nearly parallel to all the direction cosines aspossible or, what is really the same thing, a line as nearlyperpendicular to all the direction cosines as possible. In eithercase, the condition we want to satisfy as best we can is:
A * C1 + B * C2 + C3 = 0In this equation, the direction numbers of the desired plane, orits pole, are A, B, and 1. Since the equation sums to zero, welose no generality by assuming one of the direction numbers to beequal to 1 and it simplifies the math considerably.
In reality, the expression above will not sum exactly tozero but to some small value E. We want the sum of the squares ofall the values of E to be a minimum. To evaluate the sum we willneed the sums of C1*C1, C1*C2, and so on for all the directionsin our data set. Let us define a small array S such that S(1,2)equals the sum of C1 * C2, and so on. We can write a shortroutine as follows:
10 DIM S(3,3)
1000 INPUT C1,C2,C3
1010 IF C1>1 THEN 2000
1020 REM DUMMY VARIABLE TO QUIT INPUT
1030 REM COULD INPUT SOME OTHER WAY AND
1040 REM COMPUTE DIRECTION COSINES
1050 S(1,1) = S(1,1) + C1 * C1
1060 S(1,2) = S(1,2) + C1 * C2
1070 S(1,3) = S(1,3) + C1 * C3
1080 S(2,2) = S(2,2) + C2 * C2
1090 S(2,3) = S(2,3) + C2 * C3
1100 S(3,3) = S(3,3) + C3 * C3
1110 GO TO 1000
2000 REM END EXAMPLE
In this example, note the use of a dummy variable in line1010 to exit the input routine. The dummy variable was chosenbecause no real direction cosine would ever be greater than 1,and the variable number of data inputs precludes a loop of fixedlength.
The condition we need to attain is the following:
SUM OF E*E =
A*A*S(1,1) + B*B*S(2,2) + S(3,3) + 2*A*B*S(1,2) +
2*A*S(1,3) + 2*B*S(2,3) = MIMIMUMIn this case, we don't know A and B, so we treat them as thevariables and differentiate the above formula with respect to Aand B. We obtain the two conditions:
A * S(1,1) + B * S(1,2) + S(1,3) = 0
A * S(1,2) + B * S(2,2) + S(2,3) = 0We can then solve these two equations for A and B. The followingroutine will do the job.
2010 REM SOLVE FOR CONSTANTS
2020 A = S(1,2)*S(2,3) - S(1,3)*S(2,2)
2030 A = A/(S(1,1)*S(2,2) - S(1,2)*S(1,2))
2040 B = S(1,2)*S(1,3) - S(2,3)*S(1,1)
2050 B = B/(S(1,1)*S(2,2) - S(1,2)*S(1,2))
2060 C = 1/SQR(1 + A*A + B*B)
2070 A = A*C: B = B*C
2100 REM END EXAMPLEWe could have solved for the third direction cosine, C, when wederived the least-squares fit, but then we would have had a messyjob of solving three simultaneous equations. This approach is alot simpler.7-6 BEST FIT OF A CONE TO A SET OF DIRECTIONS
A circular cone is a useful approximation to a variety ofgeologic phenomena. Many folds are conical rather thancylindrical. Shatter cones, or fracture surfaces probablyproduced by shock waves during meteor impacts, have surfacestriations which define the surface of a cone. Conicaldistributions, when plotted on a stereonet, lie on a smallcircle. An application to plate tectonics might be finding thebest fit pole of rotation between two plates given a transformfault (which is ideally an arc of a small circle). Assume thecone axis has direction cosines A1, A2, and A3 and a given linewe wish to fit to the cone has direction cosines C1, C2, and C3.The angle between the cone axis and any line on the surface ofthe cone is constant and given by:
A1 * C1 + A2 * C2 + A3 * C3 = Kwhere K is the cosine of half the apical angle of the cone. Wecan simplify the equation by letting A = A1/a3, B = A2/a3, andC = -K/a3, and the equation becomes
A * C1 + B * C2 + C3 - C = 0As before, the formula will generally not be exactly equal tozero but will have a small value E, and we find the best fit byreducing the sum of the squares of all E to a minimum. We willneed the sums of C1*C1, C1*C2, and so on. The routine we used inthe previous section (Lines 1000-2000) will work just as wellhere, but with modifications because we need additional sums.
10 DIM S(3,3)
1000 INPUT C1,C2,C3
1010 IF C1>1 THEN 2000
1020 REM DUMMY VARIABLE TO QUIT INPUT
1030 REM COULD INPUT SOME OTHER WAY AND
1040 REM COMPUTE DIRECTION COSINES
1050 S(1,1) = S(1,1) + C1 * C1
1060 S(1,2) = S(1,2) + C1 * C2
1070 S(1,3) = S(1,3) + C1 * C3
1080 S(2,2) = S(2,2) + C2 * C2
1090 S(2,3) = S(2,3) + C2 * C3
1100 S(3,3) = S(3,3) + C3 * C3
1110 S1 = S1 + C1
1120 S2 = S2 + C2
1130 S3 = S3 + C3
1140 N = N + 1
1150 GO TO 1000
2000 REM END EXAMPLEWe follow basically the same approach as before. However, we havethree independent variables and the math is correspondingly morecomplex. We want to minimize
SUM OF E*E =
SUM OF (A*C1 + B*C2 + C3 + C) * (A*C1 + B*C2 + C3 + C) =
A*A*S(1,1) + 2*A*B*S(1,2) + B*B*S(2,2) +
S(3,3) 2*A*C*S1
+ 2*B*C*S2 +
N*C*C + 2*C*S3
+ 2*B*S(2,3)Again, we differentiate with respect to each coefficient A, B andC, and obtain
A * S(1,1) + B * S(1,2) + C * S1 + S(1,3) = 0
A * S(1,2) + B * S(2,2) + C * S2 + S(2,3) = 0
A * S1 + B * S2 + C * N + S3 = 0We now have three equations in three unknowns, and the easiestway to program the solution is through determinants:
2010 REM SOLVE FOR COEFFICIENTS
2020 D = S(1,1) * (S(2,2)*N - S2*S2)
2030 D = D + S(1,2) * (S1*S2 - N*S(1,2))
2040 D = D + S1 * (S(1,2)*N - S1*S2)
2050 D1 = -S(1,3) * (S(2,2)*N - S2*S2)
2060 D1 = D1 - S(2,3) * (S1*S2 - N*S(1,2))
2070 D1 = D1 - S3 * (S(1,2)*N - S1*S2)
2080 D2 = S(1,1) * (-S(2,3)*N + S2*S3)
2090 D2 = D2 + S(1,2) * (S(1,3)*N - S1*S3)
2100 D2 = D2 + S1 * (S(2,3)*S1 - S2*S(1,3))
2110 D3 = S(1,1) * (S(2,3)*S2 - S(2,2)*S3)
2120 D3 = D3 + S(1,2) * (S(1,2)*S3 -S(1,3)*S2)
2130 D3 = D3 + S1 * (S(2,2)*S(1,3)-S(1,2)*S(2,3))
2140 A = D1/d: B = D2/d: C = D3/d
2150 A3 = 1/SQR(1 + A*A + B*B)
2160 A1 = A*A3: A2 = B*A3: K = -C*A3
Return to Mathematical Formulas Index
Return to Crustal Movements Syllabus
Return to Techniques Manual Index
Return to Professor Dutch's Home Page
Created 22 January 1999, Last Update 22 January 1999